#include "core/cell-dyn.hpp"
#include "yocto/code/unroll.hpp"

#if defined(FASTMD_CELL_OPENMP)
#include <omp.h>
#endif

namespace FastMD
{
	////////////////////////////////////////////////////////////////////////////
	//
	// CellDynamics
	//
	////////////////////////////////////////////////////////////////////////////
	CellDynamicsType:: ~CellDynamicsType() throw()
	{
		
	}
	
	
	static inline size_t __get_num_threads()
	{
#if defined(FASTMD_CELL_OPENMP)
		return omp_get_max_threads();
#else
		return 1;
#endif
	}
	
	CellDynamicsType:: CellDynamicsType( const Frame &frame, const Real CutOut ) :
	CellDynamicsBase(frame,CutOut,__get_num_threads())
	{
	}
	
	
	void CellDynamicsType:: eval( Frame &frame , Potential &potential) throw()
	{
		assert( cells != NULL );
		Cells        &C        = *cells;
		const size_t  nc       =  C.count;
		const size_t  n        =  frame.particles;
		const Vector &position =  frame[ Frame::POSITIONS     ];
		Vector       &force    =  frame[ Frame::ACCELERATIONS ];
		PBC_Proc      pbc_f    =  frame.pbc_f;
		const V3D    &box      =  frame.box;
		const V3D    &ibox     =  frame.ibox;
		
		//======================================================================
		//
		// dispatch all positions
		//
		//======================================================================
		C.dispatch( node,nodes,frame[ Frame::POSITIONS ]);
		
		//======================================================================
		//
		// initialize forces/accelerations
		//
		//======================================================================
#define FLDZ(i) force[i].ldz()
		YOCTO_LOOP_FUNC_(n,FLDZ,0);
		
		//======================================================================
		//
		// process each cell
		//
		//======================================================================
#if defined(FASTMD_CELL_OPENMP)
#pragma omp parallel default(shared)
#endif
		{
#if defined(FASTMD_CELL_OPENMP)
			const size_t iSlice = omp_get_thread_num();
			assert( iSlice < cells->slices );
			const size_t c_ini = C.slice[iSlice].ini;
			const size_t c_end = C.slice[iSlice].end;
#else
#define c_ini 0
#define c_end nc
#endif
			for( size_t ic=c_ini; ic < c_end; ++ic )
			{
				Cell         &c    = C[ic];
				Cell        **link = c.link;
				const size_t  nl   = c.links;
				
				//--------------------------------------------------------------
				// initialize local potential
				//--------------------------------------------------------------
				Real V = 0;
				V3D  f_ij;
				
				//--------------------------------------------------------------
				// outer loop on inside particles
				//--------------------------------------------------------------
				CellNode *in = c.inside;
				while( in )
				{
					const size_t i   = in->index;
					const V3D   &r_i = position[i];
					V3D         &f_i = force[i];
					
					//----------------------------------------------------------
					// First inner loop: inside the same cell
					//----------------------------------------------------------
					CellNode *jn = in->next;
					while( jn )
					{
						const size_t j   = jn->index; assert(i!=j);
						const V3D   &r_j = position[j];
						V3D          r_ij(r_i,r_j);
						V3D         &f_j = force[j];
						
						pbc_f(r_ij,box,ibox);
						const Real V_ij = potential(f_ij,r_ij);
						V += (V_ij + V_ij);
						f_i.sub( f_ij );
						f_j.add( f_ij );
						jn = jn->next;
					}
					
					//----------------------------------------------------------
					// Second inner loop: over the neighbors
					//----------------------------------------------------------
					for( size_t il=0; il < nl; ++il )
					{
						assert( link[il] != NULL );
						jn = link[il]->inside;					
						while( jn )
						{
							//--------------------------------------------------
							// update V and only f_i
							//--------------------------------------------------
							const size_t j = jn->index;
							const V3D &r_j = position[j];
							V3D        r_ij(r_i,r_j);
							
							pbc_f(r_ij,box,ibox);
							V += potential(f_ij,r_ij);
							f_i.sub( f_ij );
							jn = jn->next;
						}
					}
					
					
					in = in->next;
				}
				c.V = V;
			}
		}
		
		//======================================================================
		//
		// collect data
		//
		//======================================================================
		frame.V = 0;
		for( size_t ic=0; ic < nc; ++ic )
		{
			frame.V += C[ic].V;
		}
		frame.V *= __HALF;
		
		frame.forces_to_accelerations();
		
	}
}
